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INTRODUCTION 


\* 

Atmospheric  drag  has  a  significant  effect  on  low-altitude  trajectories.  To  measure  this  force,  an  accel¬ 
erometer  will  be  tested  which  senses  this  nonconservative  acceleration.  The  information  obtained  will  then 
be  used  to  improve  the  orbit  and  obtain  a  more  accurate  gravity  field  model. 

When  the  raw  data  is  first  received,  all  known  systematic  effects  will  be  removed.  These  include  grav¬ 
ity  gradient  and  pitch-rate  induced  accelerations.  However,  the  preprocessed  data  will  still  contain  anomalies 
of  which  there  is  little  or  no  a  priori  knowledge.  Satellite  attitude  thruster-induced  vibrations  are  a  major 
source  of  unmodeled  noise.  As  for  the  instrument  itself,  each  axis  of  the  accelerometer  adds  an  undetermined 
bias  to  the  output  accelerations.  In  addition,  the  conversion  of  the  digital  output  to  engineering  units  in¬ 
troduces  a  scale  factor  error.  Other  effects  may  take  the  form  of  additive  correlated  noise  and  white  noise. 

The  first  instrument  to  be  flown  is  the  miniature  electrostatic  accelerometer  (MESA)i>uilt  by  Bell 
Aerospace  Textron  under  contract  to  the  Air  Force  Geophysics  Laboratoryif  he  MESA  senses  accelerations 
along  the  vehicle  pitch,  roll,  and  yaw  axes.  Relying  on  the  accuracy  of  the  vehicle’s  attitude  control  system, 
these  axes  are  assumed  to  be  along  the  K,  R  X  V,  and  -R  vectors,  respectively.  Hence,  any  attitude  errors 
will  be  reflected  in  the  measurements.  The  MESA  also  introduces  undetermined  scale  factor  errors  and 
biases. 

Consequently,  a  computer  program  was  designed  which  applies  a  filter  to  the  processed  accelerations. 
This  filter  should  statistically  remove  most  of  the  unmodeled  corruptions.  The  filtered  data  is  then  input 
to  a  special  version  of  the  Celest  orbit  determination  programiMn  order  to  check  out  the  filter  and  perform 
special  case  studies  prior  to  launch,  simulated  processed  data  had  to  be  created.  This  required  a  series  of 
computer  programs  whfcji  will  be  referred  to  as  the^Data  Generator."*The  individual  programs  include 
‘'Drag  Accelerations,”  lrSHost  Vehicle  Vibrations, ^“Mcrge,”  and  “Gauss-Markov.” 
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COMPUTER  PROGRAM  DESCRIPTIONS  AND  FORMULATIONS 


DRAG  ACCELERATIONS 

Updates  were  inserted  in  the  “Orbgen”  and  “Atdrag”  subroutines  of  Celest1  Version  16.  During 
integration,  the  atmospheric  drag  accelerations  are  converted  from  the  inertial  frame  to  a  perpendicular 
local  frame  (Uv  Vv  U^): 

Ui  =  RX(V>(R) 

U2  =  VXR 
U3  =  -R 

A  2-sec  step  size  is  used  to  approximate  the  MESA  data  averaging  interval. 

In  the  integration  start  routine,  there  are  multiple  iterations  of  the  first  10  time  lines.  Any  duplicate 
set  of  time  lines  must  be  omitted.  Hence,  only  the  final  set  of  accelerations  are  retained. 

The  output  file  of  drag  accelerations  is  described  in  the  section  of  this  report  entitled  “Computer 
Program  Utilization”  (p.  15).  It  can  be  used  as  input  to  either  the  Host  Vehicle  Vibrations  Program  or  the 
Gauss-Markov  Program. 


HOST  VEHICLE  VIBRATIONS 

Run  once  for  each  of  the  accelerometer  coordinate  axes,  this  program  simulates  the  response  of  the 
MESA  to  the  attitude  thruster  firings  of  a  satellite  in  orbit.  The  attitude  thrusters  are  assumed  to  be  coupled. 
This  means  that,  when  a  pair  of  thrusters  fire,  they  create  a  torque  about  the  axis  perpendicular  to  the 
thruster  plane.  The  resultant  force  is  zero,  causing  purely  rotational  vehicle  motion. 

The  magnitude  and  shape  of  the  initial  impulse  can  be  specified,  with  a  maximum  duration  of  0.5  sec. 
All  impulses  are  assumed  to  be  positive  along  the  respective  axes.  For  the  purposes  of  the  simulations  to  be 
performed,  the  direction  of  the  initial  impulses  is  inconsequential. 

Eleven  vehicle  resonances  may  be  activated  by  each  thruster  pair.  A  different  amplitude  and  damping 
constant  is  input  for  each  frequency.  Because  of  the  resolution  of  the  MESA,  a  resonance  is  assumed  to  be 
totally  damped  once  its  amplitude  is  below  10~9  g  (one  g  is  9.80621  X  10~3  km/sec2). 

The  satellite  is  assumed  to  experience  maximum  drag  acceleration  at  perigee,  the  time  of  closest  ap¬ 
proach  to  the  earth.  The  thrusters  fire  most  often  during  this  period  because  of  the  high  aerodynamic 
torques  resulting  from  the  increased  magnitude  and  gradient  of  the  gravity  force.  The  opposite  is  true  at 
apogee,  the  satellite’s  farthest  distance  from  the  earth.  At  this  point,  the  drag  acceleration  and  the  frequency 
of  thrusts  are  assumed  to  be  at  a  minimum. 

Consequently,  for  each  thruster  pair,  the  probability  of  a  thrust  is  linearly  interpolated  with  respect 
to  the  true  drag  acceleration  along  the  velocity  vector.  The  time  and  drag  accelerations  at  apogee  and  peri¬ 
gee  are  obtained  from  the  Drag  Accelerations  Program.  To  vary  the  frequency  of  the  attitude  thrusts,  the 
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number  of  seconds  between  thrusts  at  apogee  and  perigee  is  input.  The  execution  time  of  the  program  can 
be  reduced  by  increasing  the  time  between  probability  recalculations. 

To  test  for  a  thruster  firing,  a  random  number  (AO  from  the  interval  [0, 1)  is  compared  to  the  proba¬ 
bility  (P)  derived  for  that  time.  lfO<N<P,  the  thruster  is  assumed  to  have  fired. 


Frequency  Response 


Since  the  MESA  sampling  interval  is  approximately  2  sec,  the  highest  frequency  that  can  be  unambig¬ 
uously  filtered  from  the  data  is  0.25  Hz  (the  Nyquist  frequency).  An  electronic  filter  was  incorporated  in 
the  MESA  to  attenuate  any  higher  frequencies.  However,  virtual  blockage  of  the  signal  does  not  occur  until 
the  10-Hz  level.  Any  resonances  above  0.25  Hz  will  alias  into  those  below  0.25  Hz  as  a  result  of  the  size  of 
the  sample  interval. 

Simulated  accelerometer  data  must  reflect  this  condition.  Hence,  an  iteration  interval  of  0.05  sec 
(1/2 f,  wher e/=  10)  is  used  in  the  Host  Vehicle  (HV)  Vibrations  Program  to  include  resonance  at  the  higher 
frequencies.  Letting  A'  be  the  number  of  0.05-sec  time  lines  since  a  given  thrust,  the  MESA  response  at  a 
particular  vehicle  resonance  (referenced  by  ‘2’)  is 

(ATT’l  =  At  e~“ 2 ^(sin  2nftD  U(K ) 

+  2e~aiT(cos2nf^)rs[(K-i)T) 

-  e~°iT rtl(K~  2)  T\ 


where 

=  the  maximum  amplitude 
ag  =  the  damping  constant 
T  =  0.05  sec 

fs  =  the  resonance  frequency 

U  =  an  array  of  10  terms  which  determines  the  size  and  shape  of  the  initial  thruster  impulse. 

The  total  response  at  KT  seconds  after  a  thruster  fires  is  then 

Nf 

r[KT]  =  £ 

8=1 

where  is  the  number  of  resonance  frequencies  excited.  For  each  consecutive  0.05-sec  time  line,  the  re¬ 
sponse  to  the  thrust  is  computed  until  it  is  damped  below  10~9  g.  A  different  sequence  of  accelerations  is 
computed  for  each  of  the  four  pairs  of  thrusters. 
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Using  the  method  described  earlier,  the  program  checks  for  thrusts  at  specified  intervals.  If  a  thrust 
is  indicated,  the  appropriate  values  are  added  to  the  subsequent  accelerations.  The  aggregate  responses  are 
then  averaged  over  a  2-sec  period  and  time-tagged.  A  flowchart  of  the  program  logic  is  presented  in  Figure  1 
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Figure  1 .  Host  Vehicle  Vibrations  Flowchart 
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Figure  1.  Host  Vehicle  Vibrations  Flowchart  (Continued) 
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Figure  1.  Host  Vehicle  Vibrations  Flowchart  (Continued) 
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Derivation  of  Damping  Constants 

The  host  vehicle  behaves  like  a  damped  oscillator  when  excited  by  an  attitude  thrust.  The  equation 
of  motion  for  this  disturbance  is 


F(t)  -mx+cx  +  kx  =  0 

where  m  and  c  are  the  coefficients  of  mass  and  damping,  and  k  is  the  spring  constant. 

One  solution  to  this  equation  is  to  let  *  =  0(r),  such  that  x,  x,  and  x  differ  only  by  constant  multiplica¬ 
tive  factors.  The  general  exponential  form  Ae9' ,  where  A  and  q  are  constants,  fills  this  requisite  perfectly. 
Substituting  in  F\t), 


F(t)  =  m(Ae #')"  +  dA^)'  +  kiAe*") 

=  Ae^'lmq2  +  cq  +  A:]  =0 
Since  eqt  never  vanishes,  this  implies  that 

mq2  +  cq  +  k  =  0 

The  general  solution  to  this  quadratic  equation  is,  of  course, 

-c  ±  v?  -4mk 

I*  - to - 

The  real  part  of  this  root  is  the  damping  term,  whereas  the  imaginary  part  is  the  oscillatory  term. 
Since  our  system  is  characterized  by  damped  oscillatory  motion,  q  must  be  a  complex  number  of  the  form 
-a  ±  u>i,  such  that  a  and  oj  are  real  numbers.  Hence, 

-c  t  yj c2  -  4km 

q  =  -a  ±  ico  =  - 

2m 

where  c2  -  4km  <  0.  Consequently,  a  =  c/2m  and  «  =  yj 4km  -  c2  j2m. 


The  solution  is  then 


x  =  A  }e<?|f  +  A^2' 

=  +  A2e^a-iu)t 


=  e-at(A/wt  +  /t2e-'w) 

Requiring  that  x(t0)  =  0  (i.e.^4  j  =A2  =  >4/2)  and  x(tQ)  =  x0  (i.e.,x0  =A ), 
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x  =  e'at  ■  {A  12)  •  (e,wf  +  e',u>t) 


=  xae 


=  x0e~at  cos  ot 


The  number  ( N )  of  cycles  completed  at  time  /  is  n  -  cctJ2n  or,  equivalently,  t  =  2rrn/w. 

The  maximum  ratio  uf  displacement  at  time  t  to  the  initial  displacement  is 

x(t)  -at  I  j  -  . 

-  =  e  since  Icosoofl  <  1 


Substituting  for  t. 


Hence, 


x(t) 


e-2nna/u>  of 


(?)• 


-2ima 


co 


■■few?) 

-  '^Akm  -  ci  to  M0\ 

4nmc/2m  \  x0  I 
\f4km  -  c2  (  x„  \ 

ln[m) 

^r-fe) 


i 

2ir 


2  nc 


4km 


The  system  is  said  to  be  critically  damped  when  \fc2-  4 km  =  0.  The  critical  damping  coefficient 
(Cc)  is  then 


C  =  4fcm 

C 


Rewriting  n. 


n  ~ 


2rr 


1  -  (clCf 


(cICJ2 


1/2 


‘(5) 
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c/Cc  is  the  percentage  of  critical  damping  inherent  in  the  system.  If  x(f)  is  set  to  1  X  1(T9  g  (minimum 
amplitude),  n  would  be  the  total  number  of  cycles  required  to  damp  out  the  oscillation.  The  damping  co¬ 
efficient  a  can  be  determined  from  n  using  the  formula 


In  summary,  the  following  two  equations  should  be  used  to  determine  the  respective  damping 
coefficients : 

1 

n  =  — 

2n 

and 


where 

Xmax  =  *he  maximum  amplitude  of  the  disturbance 
*min  =  *he  minimum  sensed  acceleration  (1  X  1CT9  g) 
c/Cc  =  the  percentage  of  critical  damping  inherent  in  the  host  vehicle, 
a  corresponds  to  A  B  in  the  HV  vibrations  formulation, 

MERGE 

This  program  adds  the  thruster-induced  vibrational  accelerations  computed  for  each  axis  to  the  cor¬ 
responding  true  drag  accelerations  from  Celest.  If  only  one  axis  is  to  be  affected,  the  input  files  for  the 
undisturbed  axes  are  such  that  all  the  accelerations  are  set  to  zero.  These  files  should  have  the  correct  time 
tags  and  format,  however  (p.  20). 

GAUSS-MARKOV  RANDOM  NOISE 

Using  a  Gauss-Markov  statistical  process,  all  three  components  of  random  noise  are  computed  simul¬ 
taneously.  The  composite  error  vector  results  from  a  random  bias  and  scale  factor,  additive  correlated  noise, 
axis  misalignment,  and  white  noise.  The  magnitudes  of  the  errors  are  assumed  to  be  small  and,  hence,  are 
retained  only  to  the  first  order. 
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The  random  noise  is  computed  at  each  time  line.  It  is  then  added  to  the  drag  accelerations  corrupted 
by  the  thruster-induced  vehicle  vibrations.  That  is, 

as(0  =  a/t)  +  e,(r)  +  e2(t)  +  e3(t)  +  e4(t) 

such  that 

as  =  the  output  acceleration 
a{  =  the  input  acceleration 

el  =  the  error  due  to  the  random  bias  and  scale  factor 
e2  =  the  time-dependent  orientation  error 
<?3  =  the  additive  correlated  noise 
e4  =  the  white  noise. 

The  models  of  the  individual  error  sources  follows. 

Bias  and  Scale  Factor 

Much  of  the  unmodeled  system  noise  may  result  from  an  unknown  bias  or  scale  factor.  Some  types 
of  instrument  error  are  of  this  form. 

ej  =  b  +  Sa(t) 

where  b  and  s  remain  constant  for  the  duration  of  the  run.  S  =  Is-  II  where  sH;  i,j-  1,3  are  input  quantities. 

r=v+v+M- 

bx  is  randomly  chosen  from  a  Gaussian  (normal)  population  of  mean  Ubx  and  variance  o^x ,  where  p/,x 
and  Obx  are  input  constants.  In  the  same  manner,  by  and  bz  are  independently  selected  on  the  basis  of  the 
input  values p*  ol  and/uj  ,ol  ,  respectively. 

y  y  t  z 

Time-Dependent  Orientation  Error 

Since  the  satellite  has  an  attitude  sensing  and  control  system,  the  axis  orientation  errors  (0)  are  assumed 
to  be  small  angles  (i.e.,  sin  0  «=  0). 

r2un)  =  *  (t„)«(t„) 

where  a(tn)  is  the  input  acceleration  vector  at  time  /  . 
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The  nonzero  elements  of  K  (tn)  are  assumed  to  be  pair-wise  uncorrelated  for  a  given  time  tn- 


*«n> 


*jrv n' 

0 


o  M2(tn)  -My(tn) 

-■ *,«„ )  0  MJtJ 

My(tn)  -MJtn) 

such  thatjW^(rn)  are  the  current  orientation  error  angles  for  the  respective  axes 
M(tn)  itself  is  a  time-dependent  vector  described  by 

*'l>  =  VJ'O 

and 

,)  -  +  «  -  R2Jl2SmWm(‘„) 

Reflecting  the  length  of  the  correlation  period, 

"Pmv  0 


Rm  = 
m 


o 

o 


0 

0 


such  that 


exp (-Af/v-  ) 


Pmk  ~ 


At  is  the  time  between  data  points  and  rm  is  the  correlation  time  for  each  axis  k. 


H ’(tj)  is  a  vector  of  standard  normal  variates  of  mean  zero  and  o  =  1. 
m x  n1 


sm  = 

m 


0 


0 

0 


where  om.  ate  the  standard  deviations  (in  radians)  of  the  orientation  angles  for  the  individual  axes. 


Additive  Correlated  Noise 


One  known  source  of  additive  correlated  noise  is  the  effect  of  temperature  variation  on  the  instru¬ 
ment.  This  type  of  noise  is  modeled  as  the  sum  of  a  term  (fttn))  of  on-axis  correlated  noise  and  a  second 
term  ($(/„))  of  cross-axis  correlated  noise.  Both  are  obtained  using  Gauss-Markov  vector  sequences:  i.e., 

*>*>  =  Z'n)  + 
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For  every  tn,E^r\(tn)  •!} T(tn)  |  (the  expected  value)  is  a  diagonal  matrix. 

^('l>  = 

and 

*<W  = 

K  0  °' 

w =  0  %  0 

_°  0  p^_ 

where  Pflt  =  exp(-A//rt)fc).  At,  again,  is  the  time  between  data  points  and  Tr)k  is  the  correlation  time. 


aVk  is  the  standard  deviation  of  the  correlated  noise.  W  (/„)  is  a  vector  of  standard  normal  variates  of  mean 
zero  and  a  -  1 . 

For  the  second  term,  £ |F(f„)  (  =  PE>  where  is  a  nondiagonal  symmetric  input  covariance 

matrix. 

T(f,)  =  C^(f,) 

and 

T('„+1)  =  +  wv 

where 


P(k  =  exp(-Af/rtfc)  where  Af  is  the  time  between  data  points  and  r(k  is  the  correlation  time  for  each  axis. 
In  terms  of  the  input  matrix 


13 


and 


=  (Pt  -  R^)m 

,  of  course,  is  a  random  vector  like  that  of  W^. 


Additive  White  Noise 


The  uncorrelated  random  noise  is  modeled  as 


.(O 


where 


S 


W 


0 


0 


0 

0 


z 


<rWfc  is  the  standard  deviation,  or  a  priori  magnitude,  of  the  white  noise  for  each  axis.  The  elements  of  Wfn) 
are  random  numbers  from  a  set  of  standard  normal  variates  of  mean  zero,  and  a  standard  deviation  of  one. 


COMPUTER  PROGRAM  UTILIZATION 


DRAG  ACCELERATIONS  PROGRAM 

Input/Output 

The  orbit  parameters  are  left  to  the  user’s  discretion.  Reference  1  describes  the  necessary  input  values. 
However,  an  integration  step  size  of  2  sec  is  required. 

The  accelerations,  in  the  orthogonal  local  frame,  are  written  on  TAPE51.  Table  1  shows  the  format 
for  each  record. 

Table  1 .  Drag  Accelerations  Output  File  Format 


Word 

Description 

Unit 

Format 

1 

Time  from  Epoch 

sec 

Floating  Point 

2 

Drag  Along  U\  Axis 

km/sec2 

Floating  Point 

3 

Drag  Along  U2  Axis 

km/sec2 

Floating  Point 

4 

Drag  Along  U2  Axis 

km/ sec2 

Floating  Point 

The  times  of  apogee  and  perigee  are  printed  out,  along  with  the  corresponding  drag  accelerations  along 
U, .  This  information  will  be  needed  as  input  to  the  HV  Vibrations  Program. 


HOST  VEHICLE  VIBRATIONS  PROGRAM 


Input 


A  file  (TAPE2)  of  uncorrupted  drag  accelerations  must  be  input  to  compute  the  probability  of 
thruster  firings  during  the  orbit.  This  is  the  same  file  as  that  output  by  the  Drag  Accelerations  Program 
previously  described. 

The  data  cards  are  grouped  in  sections,  which  can  be  arranged  in  any  order  in  the  deck.  All  card 
identifiers  are  left-justified,  while  the  data  words  are  right-justified  in  the  field.  Table  2  describes  the  card 
format  and  Table  3  describes  input  parameters. 
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Table  3.  HV  Vibrations  Card  Input 


resonance  on  correspond-  resonance  on  correspond-  resonance  on  correspond-  resonance  on  correspond¬ 
ing  “FREQUENCY”  ing  “FREQUENCY”  ing  “FREQUENCY”  ing  “FREQUENCY” 

section  card  section  card  section  card  section  card 


Output 

The  output  file  (TAPE3)  is  a  time  history  of  the  thruster-induced  vibrational  accelerations.  Table  4 
shows  the  format  for  each  record  of  the  file. 

Table  4.  HV  Vibrations  Output  File  Format 


Word 

Description 

Units 

Format 

1 

Time  from  Epoch 

sec 

Floating  Point 

2 

Vibrational  Acceleration  | 

km/ sec2 

Floating  Point 

The  printed  output  is  arranged  in  four  sections.  First  of  all,  the  card  input  is  printed  out  and  labeled 
accordingly.  The  derived  probability  for  a  thrust  is  then  printed  at  the  specified  interval  for  each  thruster 
pair.  The  vehicle  responses  to  each  thruster  are  printed  out  at  a  0.05-sec  step  size.  Finally,  the  2-sec  aver¬ 
ages  of  the  aggregate  vehicle  responses  over  the  orbit  span  are  given.  In  addition,  the  number  of  thrusts  per 
pair  is  given  for  the  previous  2-sec  period. 

MERGE  PROGRAM 


Input 


The  input  data  card  specifies  the  epoch  of  the  accelerations.  The  format  is  shown  in  Table  5. 

Table  5.  Merge  Input  Card  Format 


Word 

Description 

Columns 

Format 

i 

Year 

1-10 

FI  0.0 

2 

Day 

11-20 

FI  0.0 

3 

Second 

21-30 

FI  0.0 

Four  binary  files  are  input  to  this  program,  fhe  first  file  (TAPE4)  contains  the  true  drag  in  the  per¬ 
pendicular  local  frame.  This  file  is  the  same  as  that  output  by  the  “Drag  Accelerations”  version  of  Celest. 

The  remaining  three  files  (TAPE1 ,  TAPE2,  and  TAPE3)  are  produced  by  separate  runs  of  the  HV 
Vibrations  Program.  For  each  MESA  axis,  there  is  a  file  of  thruster-induced  vehicle  vibrations.  These  files, 
of  course,  have  the  same  format  as  that  output  by  HV  vibrations. 


Output 

The  time  and  corresponding  corrupted  accelerations  (km/sec2)  are  printed  out  at  each  2-sec  interval. 
If  the  corresponding  times  on  two  files  are  not  equivalent,  an  error  message  is  printed  out  and  the  program 
is  terminated.  The  same  action  is  taken  when  an  end-of-file  is  reached  on  one  of  the  input  files. 
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The  output  file  (TAPES)  contains  the  true  drag  corrupted  by  the  host  vehicle  vibrations.  The  accel¬ 
erations  are  in  the  perpendicular  local  frame  in  units  of  km/sec2. 

The  first  record  has  three  words  describing  the  epoch  (year,  day,  second)  of  the  accelerations.  All  time 
tags  on  the  file  are  referenced  to  this  epoch. 

Each  subsequent  record  has  the  format  described  in  Table  6.  The  first  word  is  an  integer  specifying 
the  number  (AO  of  pertinent  observations  in  that  record.  Although  there  are  always  1000  observations  per 
record,  only  N  of  them  are  valid.  Each  observation  consists  of  four  words:  a  time  tag,  and  the  x-,  y-.  amd 
z-axis  accelerations. 


Table  6.  Merge  Output  File  Format 


Word 

Value 

Units 

Type 

r 

N 

~ 

Integer 

2 

Time  of  Observation  | 

Seconds 

Floating  Point 

3 

4 

Ax  ' 

Av  1 

H 

km/sec2 

km/sec2 

Floating  Point 
Floating  Point 

5 

^  J 

km/sec2 

Floating  Point 

6 

Time  of  Observation 

Seconds 

Floating  Point 

7 

8 

AX  1 

Av  I 

km/sec‘ 

km/sec2 

Floating  Point 
Floating  Point 

9 

J 

km/'sec2 

Floating  Point 

3998 

Time  of  Observation  ] 

Seconds 

Floating  Point 

3999 

Ax  ' 

►1000  1 

km/sec2 

Floatingpoint 

4000 

AV 

km/sec2 

Floatingpoint 

4001 

4  J 

km/ sec2 

Floating  Point 

GAUSS-MARKOV  PROGRAM 


Input 

The  input  file  (TAPE1 )  is  the  same  as  that  output  by  the  Merge  Program.  The  format  was  specified 
earlier. 

Every  input  card  must  be  labeled  with  the  appropriate  heading,  left-justified  in  Column  1 .  For  all 
17  cards,  the  format  is  (A  10,  3E20. 14).  They  are  described  in  Table  7. 
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Table  7.  Gauss-Markov  Card  Input 


Section 

Card  Label 

Columns 

Item 

Units 

General 

INTINT 

11-30 

A/ 

sec 

TIMELIM 

11-30 

Time  at  which 
to  terminate  program 
execution 

sec 

Tl 

BIASX 

11-30 

31-50 

t 

!• 

BIASY 

11-30 

31-50 

Hy 

°h 

g 

g 

BIASZ 

11-30 

Vi, 

g 

31-50 

g 

SCAFACl 

11-30 

sa.i) 

g 

31-50 

S(l,2) 

g 

51-70 

S(l,3) 

g 

SCAFAC2 

11-30 

S(2,l) 

g 

31-50 

S(2,2) 

g 

51-70 

S(2,3) 

g 

SCAFAC3 

11-30 

S(3,l) 

g 

31-50 

S(3,2) 

g 

51-70 

S(3,3) 

g 

e2 

TAUM 

11-30 

Tmx 

sec 

31-50 

Tm  v 

sec 

51-70 

Tmz 

sec 

SIGMAM 

11-30 

°mx 

— 

31-50 

°mv 

- 

51-70 

°mz 

— 

TAUETA 

11-30 

TVx 

sec 

31-50 

TVy 

sec 

51-70 

sec 

SIGMAETA 


n-30 

31-50 

51-70 


% 

a*z 


Table  7.  Gauss-Markov  Card  Input  (Continued) 


Section 

Card  Label 

Columns 

Item 

Units 

*3 

TAUXI 

11-30 

T(x 

sec 

(Continued) 

31-50 

Th 

sec 

51-70 

TG 

sec 

PX11 

11-30 

Pt(U) 

g* 

31-50 

Pt0.2) 

Pt(l,3) 

51-70 

82 

PXI2 

11-30 

Pt(2,l) 

P{(2,2) 

Pt(2,3) 

g? 

31-50 

g 

51-70 

g2 

PX13 

11-30 

Pt(3,l) 

Pf(3,2) 

g2 

31-50 

g 

51-70 

Pj(3,3) 

g2 

SIGMAW 

11-30 

g 

31-50 

0LOy 

g 

51-70 

g 

Output 


The  format  of  the  output  file  (TAPE2)  is  the  same  as  that  produced  by  the  Merge  Program.  The 
accelerations  are  written  in  units  of  km/sec2.  As  the  end  product  of  the  Data  Generator,  this  file  of  data 
is  ready  for  filtering. 

Each  data  card  is  printed  exactly  as  it  was  punched.  In  the  case  of  an  error,  an  appropriate  message 
is  given.  The  elements  of  the  time-invariant  input  and  computed  matrices  are  then  specified.  (See  Table  8.) 

Finally,  the  time  history  of  the  input  and  output  accelerations  is  given  with  the  corresponding  time 
tags.  The  data  is  specified  in  units  of  g’s  for  three  coordinates. 
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SAMPLE  CASE:  NOMINAL  CORRUPTIONS 


DRAG  ACCELERATIONS 

The  first  step  was  to  compute  a  low-altitude  inertial  trajectory  using  a  25  X  25  gravity  field.  With  one 
iteration  per  time  line,  a  tenth-order  integration  was  effected  at  a  2-sec  step  size.  The  drag  accelerations  de¬ 
rived  for  this  orbit  were  saved  for  future  use. 

The  drag  model  implemented  was  the  standard  Celest1  version.  It  may  be  noted  that  there  are  dis¬ 
continuities  in  the  plots  of  the  true  drag  accelerations.  These  anomalies  mark  the  boundaries  of  the  drag 
segments  applied. 


HV  VIBRATIONS 

A  preliminary  study  was  made  of  the  structure  of  the  host  vehicle  and  its  solar  arrays.  Assuming  a 
nonrigid  body,  six  resonance  frequencies  were  determined  to  be  below  1 0  Hz.2  These  resonances  were 
used  for  the  responses  to  the  pitch  thruster  pair.  For  simplicity,  this  was  the  only  pair  allowed  to  fire. 

The  amplitude  of  the  sensed  acceleration  for  each  resonance  was  assumed  to  be  1  pg.  To  obtain  the 
respective  MESA  output  amplitudes,  the  instrument  attenuation  for  each  frequency  was  multiplied  by  this 
value.  The  vehicle  damping  was  3  percent  of  critical.  In  Table  9,  the  actual  input  values  are  listed.  The 
thrust  impulse  itself  was  assumed  to  be  instantaneous.  Therefore,  U(  1)  =  l  was  the  only  nonzero  element 
ofU. 


The  time  interval  between  thrusts  at  apogee  was  30  sec,  while  that  at  perigee  was  1 0  sec.  The  thrust 
probability  was  recomputed  every  200  sec.  To  minimize  aliasing,  the  check  for  a  thrust  was  made  every 
1 .05  sec.  The  same  resultant  vibrations  were  added  to  all  three  axes. 

Table  9.  Nominal  Case  Resonances 


Frequency 

(Hz) 

Amplitude 

(g’s) 

Damping  Constant 

0.102 

0.97  X  10"6 

0.306137790  X  10~2 

0.591 

I 

0.26  X  10'6 

0.177379839  X  I0~] 

1.125 

0.54  X  10'7 

0.337651978  X  10-1 

1.787 

0.15  X  10-7 

0.536341408  X  10T1 

2.488 

0.64  X  10‘8 

0.746736107  X  J0_1 

3.847 

0.16  X  10'8 

0.1 15461970  X  10° 

GAUSS-MARKOV 


A  bias  of  0.1  Mg  was  applied  to  each  of  the  three  axes.  The  scale  factor  was  set  to  zero.  For  the  time- 
dependent  orientation  error,  the  standard  deviation  was  0.15°  for  each  axis.  The  correlation  time  was  one- 
fourth  of  a  revolution  or  1350  sec. 

No  additive  correlated  noise  was  applied.  However,  the  white  noise  level  was  5  X  I  O'9  g  on  all  axes. 


RESULTS 

As  expected,  the  out-of-plane  components  of  the  drag  were  much  smaller  than  the  along-track  com¬ 
ponent.  Hence,  the  plots  for  the  y  and  z  axes  were  magnified  to  expose  their  finer  structure.  The  nominal 
case  accelerations  were  submitted  to  a  Fourier  filter,  which  successfully  recovered  the  true  drag.  Plots  of 
the  simulated  MESA  data  are  attached  (Figures  2-9). 
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Figure  3.  Sample  Case:  Cross-Track  Drag  Accelerations 
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Figure  6. 
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Figure  8.  Sample  Case  Simulated  Data:  Cross-Track 


'ase  Simulated  Data:  Radial 


SUMMARY 


A  sample  case  was  run  simulating  data  with  nominal  corruptions.  Both  vibration-induced  accelerations 
and  Gauss-Markov  modeled  noise  were  added  to  the  drag  computed  in  Celest.  The  simulated  data  appears 
quite  realistic.  In  fact,  it  closely  resembles  the  real  data  received  from  the  first  orbiting  three-axis 
MESA.  A  Fourier  filter  removed  virtually  all  of  the  ocrruptions  from  the  true  drag  accelerations.  The  details 
of  the  runs  are  described  in  the  section  of  this  report  entitled  “Sample  Case:  Nominal  Corruptions ”(p.25). 

The  “Data  Generator”  is  quite  versatile  in  that  vibration-induced  accelerations  and  Gauss-Markov 
noise  can  be  added  independently  and,  if  desired,  exclusively  of  one  another.  One  can  add  corruptions  to 
any  combination  of  axes. 

In  the  Host  Vehicle  Vibrations  Program,  1 1  inherent  vehicle  resonances  can  be  designated.  These 
vibrational  frequencies  can  be  initiated  by  any  one  of  four  pairs  of  thrusters.  The  amplitude  and  damping 
constants  of  the  resonances  are  different  for  each  thruster.  In  addition,  an  initial  thrust  impulse  can  be 
superimposed  on  each  of  the  excited  vibrational  modes.  The  size  and  duration  of  this  impulse  can  be 
specified. 

Up  to  five  types  of  first-order  Gauss-Markov  noise  may  be  added  to  the  data.  The  noise  may  take  the 
form  cf  a  bias  or  scale-factor  induced  anomaly,  a  time-dependent  orientation  error,  on-axis  and  cross-axis 
additive  correlated  noise,  or  white  m  >•«*.  This  model  can  therefore  emulate  an  instrument  bias  on  each  axis, 
digital-to-engineering  unit  conversion  errors,  vehicle  attitude  errors,  temperature  variation  effects,  and  uni¬ 
formly  distributed  random  noise. 
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FORTRAN  LISTINGS 


DRAG  ACCELERATIONS  PROGRAM 


•iUcNT  AGitU 
•INitRI  ORo.»c.N.i 

GOMiON/ AGJ  tLS/i  GO*  I  ri,Ti^,A2i.,ac^ 

♦INitKT  ORObtN.29 
Hit- it 

prut  jita 
Jj*b  FORXAl  (HI) 

•GO.HPlLfc  ORrfjt  N 
►iNitRT  A  f  ORAn.  i 

JiMcNblOH  A I NEK I  1^,11 , APcRMitl ) 

Rt al 

oOHrt  Ji^/AoJ- LS/IGU,  i  Tl,r  T2,A;i,AG2 
•iNStKr  ahjkau.ia?’ 

»  * 

*  ONIT  JUPLtCAT:.  TINE  LINtS  * 

*  • 

IF  (IbU»c.R«i*  OO  TO 
IF  (Ti.EG.U.i  jO  (0  9U0 
IF  (Tl.cT. 2u.)  vj  o  (0  9u  9 
ItiOsIbOfi 
OU  r J  9u9 

*  * 

*  OcRitfc  MAIRIX  FOR  CUNt/ERSIQN  FKON  INEkTIAL  TO  LUGAl  rRME  * 

»  » 

9ua  KMAi  =  S».4Kr(X(ll*X(l)*X(2>*X(2»*X  (il»X(3>  ) 

rfMA  j- jjR  HUll)  "AO (l)+<0(2l»X0(2)tXD(i>  *XD(3I) 

9XRX=Xli)  *X0(2J  -X(2)*X0(3» 

I/XRT=X(1>*XJ(3J  -X(il»XJ (1> 
t/XR2  =  a  (  2  )  *  X  J  (  1)  -  X(  1)  *X3  (2) 

rfXKNAi,^  jaxr  (  *ARX*t/AKXt¥ARir*rfX<f  ♦rfXRZ’tfARZ) 

RX>MA.»=tf  XR*Ai» 

GOii>Art=(X( I>  *X0 ( 1) *X<2>  *XJ<2>  *X ( 3)  *XJ C 3 ) > / ( RMAG *VMA3> 

SINiAa  =  KXi/NAb/(RriAb*VMAbl 

&or»A.i=(;j>iiArt/biNiAH 

GSG  3  AH  =!./:>  INGAM 

AIN.RI  <i.,i»=UX 

AINERI (2* 1) -JV 

AINERI  (3. 1 1 -02 

ITQPlit 1»  =  XOU)  *G:>GGXM/»IMAi-X(i»  *C0Tv»4M/RMAG 
ITOMlf  21  =  xO(2>*GiGGAM/VNA;;-<(2)*;;or.,AM/RMAi 
ITG3(i»3)=X0ti)  *USGGAM/  t/MAG-X(3)*G0TGA  M/RMAG 
IIOP(2.i>  =  WXRX/tfXRMAb 
IT0P(2»2)  =  llXRy/il  XRMAG 
IT0P(2#i) stfXRZ/tfXRNAb 
I TOP (if  1>=-1.*X(  D/KMA* 
lTUP(it2)=-l«*X(2) SRHAG 
1  TOP  (if  i)  =  -A.*X(i)  /RMA* 
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*  » 

*  ,;oM*£Rr  inertial  accelcRAT  iqn;>  ro  local  frame  » 

*  * 
CAL.  HPROJ  (3,3,l,ITUP,3,AlN-Rr,3,APc.RP,3) 

HrtiTt  <91>  li  .  (APERP  (1,11,1=1  ,3)  «  (X  (  I)  «  1=1, 3>  ,  <XO(I),  1  =  1,  I) 

*  * 

*  DETERMINE  APJGEE  AND  PERIGEE  * 


it  cu.ur.2.1  co  ro  3353 

IF  (TI.Ea.2.1  60  TO  3331 
TT1=0.  3  ACl=APc.RP  4 1«  II  S  30  TO  909 
3331  TT2=2.  $  AC2=APERP Cl. II  S  30  TO  909 
3333  IF  (AC2.LT .ACi.ANO.AC2.LT.APER9 (1,1>»  PRINT  111,  TT2.A02 
IF  ( A02 .6T ■ ACi. ANQ.AG2. GT. APtR9 (1,111  PRINT  112,  TT2,AC2 
TT1*TT2  3  AC1=AC2 
TT2-TI  3  AC2-APERP  (1, 1) 

60  TO  9U9 

111  FORMAT  (H0.5X,*PERI6£E  —  TINE  *  *,F».  0,10X,*ACCELER1TI3<(  =  *, 
IE  20.141 

112  FORNAT  (  H 0,9X,* APOGEE  —  TINE  =  • ,Fb. 0 ,10X, ♦ACCELERATION  =  », 
1E20.1I*> 

909  CONTINUE 
•COMPILE  AT  DRAG 
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HOST  VEHICLE  VIBRATIONS  PROGRAM 


PROGRAM  H7VIBE  ( INPUT  ,  QU  TPUT ,TAPE1=INF(JT,T  APE  2*  TAF,£3»TAPE4> 
DIMENSION  MAXS(4)  ,  SUM  (•»,  .,3 , 25  1) 

0  IMENSION  AMP(4,11)  ,0AMP<4,11)  ,FREQ(11>  ,UM  (4,  U> 

DIMENSION  T(20)  ,AC(2C),ThR1(2C>  ,THR2(2u)  ,THP(2G)  ,THY(2  3),MRH22) 
DIMENSION  MR2(22) ,MP<22) 

DIMENSION  MY (22 1  ,8P 1(221 ,3S2 (2  2)  , PP (2 2 > , B Y ( 22) , T TES T (22) . ACC(3, 1) 
DIMENSION  RESP(lCOi!C)  ,M  AX (h) 

DIMENSION  PERTH  1C >  ,PERAC(iC)  .APOTHlu), APOACC10» 

REAL  NOMTLbT 

*  * 

*  READ  CARD  INPUT  * 

*  » 

oo  5  in, li 

RE  AO  (1,10)  CAROIO. CCNST1,C0RST2,C0NST3,C0NST4 
10  FORMAT  (Alii,£15.9,E15.9,tl5.9,E15.9) 

IF  (CAROIO. NE. lGHNOMdaTtST  )  Go  TO  IOC 
NCMTL8T=C0NST1  S  GO  TO  5 
ICC  IF  (CAROIO.NE.lOHAMPLITUoES)  GO  TO  131 
DO  15  J=l,ll 

15  READ  (1,13)  CAROID, <4ilP( JJJ, J) , JJJ=1, h) 

GO  TO  5 

131  IF  (CAROIO. NE.IOHOAPFING  )  GO  TO  132 

DO  lb  K=  1 » 11 

lfe  READ  (1,13)  CAROIO, (0  AMP (KKK ,K ) ,KKK  =  1,4) 

GC  TO  5 

1C2  IF  (CAROIO. NE.IOHPPOOCALINT)  GC  TO  133 
Pfi03INT=C0NSTl  $  GC  TO  5 
103  IF  (CAROID. NE.10HT8TFAPO  )  GO  TO  134 

T3TFAR1=C0NST1  f  T BTFAR2=C0NS T 2  $  TBT F AP=CON ST  3 

T0TFAY=CONST4  i  GO  TO  5 
1 C 4  IF  (CAROID. NE.10HTBTFPERI  )  GO  TO  105 

TBTPR1=C0NST1  t  T BT FPR2=CONST 2  $  TDTFPP=C0NST3 

TBTFPY=C0NST4  t  i  0  TO  5 
ICS  IF  (CAROID. NE. 1CHPEFACCELS  >  GO  TO  13b 
READ  (1,10)  CAROIO, PN 
NP=PN 

DO  70  M= 1 , N P 

7  L  RE  AO  (1,10)  CARDIO,  PERTH  Ml  »PERAC(M) 

GC  TO  5 

10b  IF  (CAROIO. NE.IOHAPOACCELS  )  GO  TO  107 
READ  (1,10)  CAROIO, AN 
N  A  =  AN 

DO  71  M=1,NA 

71  READ  (1,13)  CARDIO, APOTI (M> ,APOAC(M) 

GC  TO  5 

107  IF  (CAROID. NE.10HFREQUENCY  )  GO  TO  138 
00  2b3  JJJ=1,11 

263  READ  (1,13)  CARD  10 , FREQ (JJJ ) 

GC  TO  5 

ICtt  IF  (CARDID.NE.10HUM  ARRAY  )  GC  TO  109 
00  264  JJJ=1,1 0 

264  READ  (1,13)  CARO  10,  (UM < I II, JJJ ), 1 11  =  1 , 4) 

GO  TO  5 
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...  .J**'  • 


i  IF  (CARDIO .NE.10HTI ME  LIMIT)  GO  TO  110 
TIM£LIM=C0NST1 
GO  TO  5 

I  PRINT  3,  CAROIO 

I  FORMAT  (IX, ’ERROR  ON  DATA  CARDS - ID  =  *,A10) 

>  STOP  “ERROR  ON  DATA  CARDS*' 

»  CONTINUE 

PRINT  CARO  INPUT 

PRINT  200 

PRINT  901.  NCMTL8T 

L  FORMAT  (5X,*NO.  OF  TIME  lInES  BETWEEN  THRUSTER  CHECKS  -  *,F5.0) 
PRINT  911,  ( (I, FREQ  (I)) ,1  =  1,11) 

.  FORMAT  (///,5X, 'FREQUENCIES*,/,  «/ ,5  X,  12,  10X.E15. 9)) 

PRINT  90  3 

i  FORMAT  (///,5X,* AMPLITUDES*,//) 

OO  634  J JJ=  1,4 

t  PRINT  902,  JJJ.UMPl  JJJ,  J)  ,J=1,11) 
l  FORMAT  (1H0 ,SX,I1,3  t/,5 ( SX,£15,9))> 

PRINT  904 

,  FORMAT  (/// ,5X,*DAMF1NG  CONSTANTS*,//) 

00  635  KKK=  1 «  4 

i  PRINT  902,  KKK,(OAMP(KKK,K),K=l,ll) 

PRINT  905,  PROSINT 

>  FORMAT  <///,5X,*PRO0ASILITY  CALCULATION  INTERVAL  -  *,F6.0) 

PRINT  200 

PRINT  906 

>  FORMAT  (  /// ,5  X, *T  IME  3-TwEEN  THRUSTER  FIRINGS*) 

PRINT  9 07,  T8TFAR1, T5TFAR2 ,TC TF AP , TBT F AY 

’  FORMAT  (/,5X,*A30GEE*»/» 

15X  ,*R0LL1  -  *,F6.0,i:x,*ROLL2  -  **F6. D,1CX,*PITCH  -  *,F6.0,1CX, 
2* YAH  -  * , F  6 . 0 ) 

PRINT  90  2,  TBTFPR1 ,  TbTFFR£,Tt>TFPP»T3TFPY 
i  FORMAT  (/,5X,*PERIbES*,/, 

15  X  ,*R0  LL 1  -  *,F6.0,  lvX, ’ROLL  2  -  *,Fb. fl  ,1G X, *PITC«  -  *,F6.0,10X, 
2*YAW  -  *  ,  Fb  ,  C) 

PRINT  909 

i  FORMAT  I///, 5X, ’ACCELERATIONS  AT  PERIGEE*,//) 

DO  535:  1=1, NP 

;  PRINT  72,  PERTI(I) ,PERAC<I) 

:  FORMAT  ( 1H 0 , 4X, *T IME  =  *, F6. 0 , 1 0 X,* ACCELERATION  =  *,E15.9> 

PRINT  910 

FORMAT  <///,5X, ’ACCELERATIONS  AT  APOGEE*,//) 

00  5025  1=1, NA 

>  PRINT  72,  APOTKI)  ,  APCAC(I) 

PRINT  200 

PRINT  912 

■  FORMAT  ( /// ,5X, *UM  ARRAY*,/) 

00  6<, 0  111=1,4 

;  PRINT  902,  III,  (UM(  III,  JJJ),  JJJ  =  1,13) 

PRINT  913,  TIMELIM 

5  FORMAT  (///,5X,*TIME  LIMIT  -  *,F6.0I 
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COMPUTE  PeOBAdUITI  £S  FCR  THRUSTER  FIRINGS 


i 


*** 


PRINT  200 
Zli  FORMAT  (1H1I 

TMPT£ST=.05*NOMTL8T*.C5 

♦  ♦ 

*  JtRlYE  THRUST  FF J 9A3  11  IT  IE S  AT  APOGEE  ANO  PERIGEE  * 

*  * 

ROU.im=TMPTEST/T0T  F  =  Rl 
R0LLiX2=T1PTEST/TtiTFPR2 
ROLLMNl= TMPTEST/To  TFAR1 
ROLL 1N2= TiP TEST /TBTFAR2 
PICrMA  x=  TIP  TEST/  T8TFFF 
PICH1lN=TMPT£ST/TbTF4P 
VAw.1AX  =  TMPTEST/ToTFPY 
YAh.-1IN=TMPTEST/TBTF  AY 

k=  :• 

IF  (PERTH!)  .LT  .APOTH1M  GO  TO  7-2 

00  73  1=1, NA 

K=KH 

T (K»=AP0TI( I) 

AC (< I  =  APOAC (II 
THRl <<)=  R0LLMN1 
TPR2 (K )  =  P JL  LMN2 
T  PP ( K I *P ICHMIN 
THY(KI  =YAK'IIN 
KM AX=K 
<=<♦1 

IF  (I.GT.NPI  GO  TO  7-.  3 
KMA<=K 

T  <M  =P  :RTI  (  II 
AC (<) = PE  RAC  (  I  I 
T  PF 1 (< )=  R3  L  LM  XI 
Th^2( <I=R0LLMX2 
THPUI  =?ICrlMAX 
7  3  T  M  Y  (  <  I sYAHMAX 
GC  TO  7-3 
7-2  00  7-  1=1, N® 

<  =  K»1 

T (Kl =PERTI( II 
AC(XI=FE=AC(II 
TPkI <KI=R0LLmX1 
TH-:2(KI=R3lLMX2 
T HP ( < I =P IC  HM A  X 
TPY(K)=YAHMAX 
KM AX=K 
<=<♦1 

IF  (I.GT.MA)  GO  TO  7-3 
KM A  <  =  K 

T  (Kl  =  APQT  I  (  1 1 
AC<<M4PJAC<II 
THR1(K)=R0lLMN1 
TPR2(K  1=  RO  L  LMN2 
THP( K) =PICHMIN 
7-  ThY  (  <)  =Y  AM  M  IN 
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Ini  KMAXP1=KMAX*1 
KMAXP2=K'1AX*2 
KM  AXM 1=KM AX  - 1 
O AXH2=<MAX-2 

♦ 

*  CGMPUTE  THRUST  PROBABILITY  VS  Or,  AG  ACCC  Lc.  RAT  ION  FOk  EACH  HALF  REV 

♦ 

IF  ( T  ( 1)  .Nt.C.)  GO  TO  7.10 
00  7J00  1=1,  KMAXMi 
0=1*1 

MR1 (I)=(THR1 (J)-THRi (I) ) / (AC ( J) - AC ( I > ) 

MR2(I)=(THR2( J)  -THR2(  I) )  /  (  AC  (  J)  -  AC{  II  ) 

MP  f I »= (T  HP ( J) -THP ( I  I ) / t  AC( J ) - AC ( I ) I 
MV (I )= (THY ( J) -THY ( I ) ) / ( AC ( J) -AC (  I) ) 

BR1 ( II =fHRl ( Jl -MR1 ( I) *AC( J) 
aR2<II=THR2(JI-MP2(I)*AC(J) 

8P(I)=THP( J)-MP( I) * AC( J) 

7  U  GO  8Y  (II=THY(JI -MY (I) *  AC ( J ) 

MR1 (KMAX )  =  MR1 (KMAXM2) 

MR2(KMAX)=MR2(KMAXM2) 

MF  (KMAXI =MP (KMAXM2) 

MY (KMAXI  =MY (KMAXM2) 

3R1(KMAXI=THR1(KMAX  l-MRKKMAXI*  AC  (KMAX) 
3R2(KMAXI=THR2(KMAX)-MR2(KMAX) *AC(KMAX) 

8P( KMAXI =THP(KMAX) -MP ( KMAX) * AC ( KM AX) 

3Y  (KMAX) =THY (KMAX) -MY (KMAXI* AC (KM AX) 

MM=KMAXP1 
00  7J01  1=1, KMAX 
7  j  01  TTEST (I) =T ( I ) 

TTE3T(KMAXP1 1=999999, 

GO  TO  7350 

7011  00  702u  1=2, KMAX 
J=I-1 

MR1(II=(THR1(II  -THRKJ)  I  /(AC(I)-AC(  J)  ) 
MR2(II=(THR2(I)-THF2(J))/(AC(I)-AC(JI I 
MF(I)=(THP(I)-THP(jn/(AC(I)-AC(J)> 

7u2t  MY  < I )= (THV (I) -THY (  J ) ) /( AC( I)  -AC ( J)l 
MRl(l)  =MR1  (3) 

MR2 ( II =MR2 (31 
MF ( 1 I =MP (31 
MY  ( 1)  =  MV (31 

MfiKKMAXPl  >=MR1(  KMAXMI) 

MR  2  (KM  AX  PI  I  =MR2  (KMAXMI) 

MP  (KMA X3 1 )  =MP  (KMAXMI) 

M  V  (KM  A  XP  11  =  MY  (KMAXMI) 

00  7021  1=1, KMAX 

8R1(I) =THR1(II-MR1( I)*AC(I) 

BR2(I)=THR2(I)-MR2(  I) *AC(I) 

BP)I>=THP(I)-MP(1)*AC(II 
7  j  21  BY  (I)  =  THV(I )-MV (I)*AC(I) 

8RKKMAXP1I  =THP1(KNAXI-MR1(KMAXP1) *AC (KMAX) 

8R2 ( KM AXP1 1 =THR2 (KM  AX ) -MR2 (KMA XP1)* AC (KM A  XI 
BP  (KMAXP1)=TMP(KMAX)-MP(KMAXP1I*AC(KMAX) 

BY(KMAXPl)  =  T  )-Y(  KMAX  )  -MV  (KM  AXP 1 )  *  AC  ( KM  AX) 

MM=KMAXP2 

TTEST ( 1) =0 ,  I  TTEST (KMAXP2I  =999999. 

00  7 J22  I  =  2,KMAXP1 


J=I-1 

7  a  22  TTEST(I» =  T(J» 

7j5L  PRINT  dOl 

8 Cl  FORMAT  <13X,49H***  PROBABILITIES 
12H**,/// ) 

TlAST=-2. 


THRUSTER  FIRINGS  PER  SECONO  *, 


DETERMINE  NUMBER  OF  PROBABILITY  SEGMENTS 
KIIM=TIMELIM/PR0BINT+1 

COMPUTE  PROBABILITY  OF  RCLL»  PITCH,  AND  YAW  THRUSTER  FIRINGS 

FOR  EACH  INTERVAL 

00  12  K=l*  KLIM 

KK=K*1 

RK=K 

DETERMINE  THE  CENTER  TIME  LINE  CF  THE  INTERVAL 

RCH=(PR08INT* (RK-.5  >)/2 . 

RCK=AINT (RCHI 

IF  ((RCH-RCKI.GT..5)  GO  TO  163 
TIflE=2.*RCK  $  GO  TO  17 
163  TIME=2.*RCK*2. 

17  CONTINUE 

ISKZP* (TIME -TLAST) / 2-1 
IF  (ISKIP.EQ.OI  GO  TO  6343 


» 

* 

* 

* 

* 

* 


* 


4  ' 

* 

¥ 

?  •  1 

f 

* 

* 

SKIP  TO  THE  APPROPRIATE  TIME  LINE 

¥ 

¥ 

j 

JO  6342  J=1 * ISKIP 

READ  (21 

IF  ( EOF ( 2) )  25,6342 

i  * 

63<«2 

CONTINUE 

¥ 

f  1 

i 

♦ 

READ  THE  X-AXIS  ACCELERATION 

¥ 

¥ 

63<*3 

RE  AO  (2)  TI,(ACC(I,1>,I=1,3> 

IF  (E0F(2»)  25,150 

m 

15C 

IF  (TIME.EQ.TH  GO  TO  675 

¥ 

♦ 

PRINT  AN  ERROR  MESSAGE  IF  THE  TIMES  OONT 

MATCH 

• 

¥ 

PRINT  300,  K, TIME, TI 

3  JC 

FORMAT  (IX, *K  =  * , I 5, 1H, , 10 X ,*T IFE  *  *,F6.0,1H,, 

10X,*TI  = 

•*F 6*01 

* 

2666 

STOP  2666 

6  75 

ACCEL=ACC(1,1> 

TLAST=TIME 

¥ 

¥ 

COMPUTE  THE  PROBABILITY  OF  ROLL,  PITCH,  AND  YAW 

THRUSTS 

¥ 

♦ 

4^ 

FOR  THE  INTERVAL 

¥ 

¥ 

00  47  1*1, MM 

II=I*i 

IF  (TTESTm.GT. TIME. OR. TIME. GT.TTESTIIIM  GO  TO 

47 

PR08R1=BR1(I»*MR<I»*ACCEL 
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PR0BR2=3R2 ( I) ♦HR2  1 1 ) * ACCEL 

PR08P=BP ( I ) *MP( I > *  ACCEL 

PR08Y=BY (I ) ♦MYl  I ) * ACCEL 

WRITE  14)  PROBR1 , PRQ8R2 « PROBP .PRO BY 

PRBR1=PR03R1/TMPTEST 

PR8R2=  PR0BR2/ TUP  TEST 

PR  BP=PRO  BP  /  TMPT  £  S  T 

PRBY=PR08Y/TMPTEST 

PRINT  800,  K,TIME,FRBR1,PRUR2,PR8P,PRBY 
322  FORMAT  ( IX  » I 4,5X  ,*TIM£  =  * ,F E. 0 , 1H, ,5X,*R OLL1  =  *,F8.5,1H, ,5X, 
2*ROLL2  =  *,F8.5,1H,  ,5X, 

2*PITCH  =  *,F8.5,1H,  ,5X,*Y AW  =  *,F8.5I 
GC  TQ  12 
<♦7  CONTINUE 
5668  STOP 
12  CONTINUE 
GC  TO  27 

26  CONTINUE 

27  REWIND  2 
REWIND  4 


PRINT  20 G 
PRINT  803 

803  FORMAT  l 15 X , 24H***THRUST£R  RESPONSES***.///! 

PI=AC0S<-1. ) 

DO  21  1=1,4 
21  MAXS(I)=0 
DO  25  1=1,4 
DO  23  J=  1,  40 
00  23  K=  1 »  2 5 1 
23  SUM(I,J,K!=0. 

•  * 

*  COMPUTE  THE  MESA  RESPONSES  TO  THE  INDIVIDUAL  THRUSTER  INPULSES  * 

*  * 

DO  1U1  1  =  1  ,4 
DO  514  JJ=1, 10100 
514  RESP(JJI=0. 

MXOL  0=0 

*  * 

*  COMPUTE  THE  RESPONSES  AT  EACH  RESONANCE  FREQUENCY  * 

*  * 

DO  1122  J=l,ll 
RLMK1=G.  S  RLMK2=C. 

00  1133  K=l, 10003 
IF  (K.GT.10I  GO  TO  45 

E1*AMPCI,JI*EXP(-,05*OAMF(I, Jl I *SIN ( .  i*PI*FREQIJ) )*UH(I,K) 

45  £2=2. *EXP(-. 05*0 AHP  f I,J» l*CO S(. i*PI*FR£Q< J) )*RLMK1 
IF  (FREQ(J) .NE.5.1  GO  TO  2C 
E2  =  3. 

20  E3=-1.*EXP (-.1*0  AMP (I,J)>  *RLMK2 
RLM=E1*E2*E3 
RLMK2=RLM<  1 
RLMK1=RLN 

RESP(K)=RESP(K)*RLM 

MX=< 

ABRLM= ABS(RLM) 


44 


TH£  RESPONSE  IS  ASSUMEC  TO  3E  TOTALLY  0  AMP  EO  8EL0W  l.XlO**-9  GS 


IF  (A8RLN.LT. i.E-09)  GO  TO  1346 
1133  CONTINUE 

PRINT  418,  RLMK1, I, J 

Lie  FORMAT  (1H1,1X,*RLM  ARRAY  TOO  SMALL - RLMK1  =  * , E15 . 9, 2X , ✓, IX 

1*THRUST£R  *,  12, 5X,  'FREQUENCY  *,121 
3666  STOP  3666 


COMPUTE  THE  NUN9ER  OF  .C5  SEC  TIME  LINES  RE3UIRE0 
TO  0 AM P  THE  RESPONSE  TO  THE  THRUST 

1346  NAX(I>=MAXO(MX,NXOLO> 

MX0L3=MAX(I> 

1122  CONTINUE 

MXW=MAX(  I) 

PRINT  859,  I, <RESP< III) ,111=1 ,PXWl 
850  FORMAT  < INI ,5 X, I 1 ,////,(/, 6(5 X , E 15 . 9>> > 

DEPENDING  UPON  THE  POSITION  OF  THE  THRUST  WITHIN  A  2  SEC  INTERVAL 
AVERAGE  THE  APPROPRIATE  RESPONSES 
FOR  EACH  SUBSEQUENT  2  SEC  TIME  LINE 

CO  63  IS=1,4Q 
ADD- 3 .  $  KKS=1 


IS  REFLECTS  THE  POSITION  OF  A  THRUST  IMPULSE 
WITHIN  A  TWO  SECOND  INTERVAL 

DO  22  JS=1,IS 
ADO=ADD*RESP(JS) 

IF  ( JS.EQ.MAX  (III  GO  TO  53 
22  CONTINUE 

SUM( I, IS , 1 ) =A9D 
JS=IS 
40  ACO=tJ. 

KKS  COUNTS  THE  NUMBER  OF  TWO  SECCNO  AVERAGES  FOR  THRUSTER  I 
JUS=0  3  KKS=KKS*1 

JS  COUNTS  THE  NUMBER  OF  RESPONSES  AOOEO  FOR  THRUSTER  I 
*»6  JS  =  JS*1 

JJS  COUNTS  THE  NUMBER  CF  RESPONSES  AuDEO 
FOR  THE  CURRENT  TWO  SECONC  INTERVAL 

JJS=JJS*1 
AQO= AOO'RESP ( JS ) 

IF  (JS.EO.MAXdn  *0  TO  50 
IF  (JJS.NE.40J  SU  TO  4b 
50  SUM  (I  ,IS, KKS»=AOD 

IF  ( JS.NE,MAX(I>  »  GO  TO  *,0 
11  MAXS(II=MAX9 (MAXS(I) *  KKS  > 

6l  CONTINUE 


1111  CONTINUE 

MAXS1=MAXS ( 1> 

MAXS2=MAXS<2) 
mAXS3=MAXS  1  3* 

MAaS4=MAXS14* 

MX=MAXC<MAXSU*,MAXSI2)  ,  HAXS  (3  I  *  FAXSI  4** 

M*X=MX-1 

*«****♦♦«♦♦#»**►***##**»****»**»»***•.**••*****♦*•**********♦*»•******»*• 
♦  * 

*  COMPUTE  THE  AGGREGATE  TMRUST-INOOCEO  ACCELERATION  * 

*  FOR  EACH  TWO-SEC  TIME  LIKE  OF  THE  TRAJECTORY  * 

*  * 

PRINT  20G 
PRINT  880 

88«  FORMAT  ( 15 X , 31H*** V IBk ATIONAL  ACCELERATIONS****///* 

NCMTL*N0MTL8T 

IfiCLLl=IR0LL2=IPITCH=IY AW*t 
OC  518  11=1, MX 
518  RE SP  ( II)  =0  . 

RNOSET  =0  • 

INT=-1 
IPT=-1 
TM=- . 0  5 
IPIFST  =-2 
44fc  INT=I NT* 1 
TH=TM*  .05 

IF  (TM.Ea.0.1  GO  TO  6uP 
IF  1InT.LT .40*  GO  TO  4756 
600  RE$P<l)=RtSP(l*/4Q. 

FRINT  87  Q ,  TM,RESP(1> 

870  FORMAT  <5X,*TIME  =  *  ,  Fo  .  C  ,  1C  X,*  ACCElER  ATION  =  *,  f  15 . 9) 

PRINT  650,  IROLL1, I ROLL2 , IRITCH , 1YAW 
650  FORMAT  I5X , ’THRUSTER  FIRINGS  1  RCLL1  =  * , 1 2 , 1 JX , *F OLL2  =  *,I2,1jX, 
l’PITCH  =  *»I2*1CX»*  YAt,  =  * ,12 , /  ) 

*  * 

*  CONVERT  ACCELERATIONS  TO  UNITS  OF  Krf/SECE  AHO  WRITE  ON  TAPE3  * 

,  * 

RE  SP<ll=RESP«l**.0T93f'621 
WRITE  (3)  T  M  »  RE  5  P  ( 1  > 

QO  3462  1=1, MXX 
J=I*1 

3462  RESP(I*=RESPU* 

RESP (MX*  =3  . 

IR0LL1=IR0LL2=IPITCH=IYAW=C 
I  NT=t5 

*  * 

*  CHECK  IF  TIME  LIMIT  HAS  SEEN  REACHED  * 

♦  * 

4756  IF  tTM.GT.T IMELIMI  GO  TO  4866 
IPT=IPT*1 

IF  tTM.EQ.Q.I  GO  TO  5C0 
IF  TIPT.NE .NOMTL*  GO  TO  446 


CHECK  FOR  THRUSTER  FIRING 


5u0  TI2=TM/2. 

TI=AINT(TI2> 
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NS  REFLECTS  THE  TIME  OF  THE  THRUST 
WITH  REFERENCE  TO  THE  OVERLYING  TWO-SEC  INTERVAL 

NS*20.*(2.-TM*TI*2.  > 

IPI=T  M/P R36INT*1 

IF  (IPI. £3 • IFIFST I  60  TO  333 

RE  AO  (4)  PR08R1, PR0BR2, PRObP.PROBY 

IP IFST=IPI 

DETERMINE  IF  THE  ROLL  THRUSTERS  HAVE  FIRED. 

IF  SO.  ACO  THE  APPROPRIATE  ACCElcR ATIONS  TO  THE  TOTAL 
AT  EACH  OF  THE  SUBSEQUENT  TrfO-SEC  TIME  LINES. 

333  RNOSET  =RNDSET *2 . 

CAcL  RANSET  (RNDSET) 

P=RANF (DUM) 

IF  (P.GT .PR0BR1I  GO  TO  Go 
00  64  1= 1 1 M AX Si 

64  RE$P(IMRESP(mSUM(l,!)S,  I) 

IRCLLl=IROLLm 

66  RNOSET =RN3SET 42. 

CALL  RANSETIRNOSETI 
P=RANF(OUM> 

IF  (P.GT . s  R0BR2)  GO  TO  65 
00  6 7  1= 1 . M AXS2 

67  RESP(I)  =  RESP(I»4SUM(2,NS,I) 

IR0LL2=IR0lL2*1 

DETERMINE  IF  THE  PITCH  THRUSTER  HAS  FIRED. 

IF  SO,  A 00  THE  APPROPRIATE  ACCELERATIONS  TO  THE  TOTAL 
AT  EACH  OF  THE  SUBSEQUENT  TWO-SEC  TIME  LINES. 

65  RNOSET  =RNOSET  *2. 

CALL  RANSET (RNOSET) 

P=RANF (DUMI 

IF  (P.GT .PROBP)  GO  TO  343 
00  342  I=1,MAXS3 

342  RESP(I)=RESP(II +  SUM  (3.NS, I) 

IPITCH=IPITCH*1 

DETERMINE  IF  THE  YAW  THRUSTER  HAS  FIRED. 

IF  SO,  ADO  THE  APPROPRIATE  ACCELERATIONS  TO  THE  TOTAL 
AT  EACH  OF  THE  SUBSEQLENT  TWO-SEC  TIME  LINES. 

343  RN0S£T=RNDSET*2. 

CALL  RANSET (RNOSET) 

P- RANF (OUM) 

IF  (P.GT.PROBY)  GO  TO  468 
00  467  I«1 , MAXS4 
h 6 7  RESP(I>=  PESP(I)4SUM  (4, NS,  I) 

IY  AM=IV AW4 1 
468  IPT=-1 

GO  TO  446 
4666  STOP  4666 
ENO 


MERGE  PROGRAM 


PROGRAM  MERGE < INPUT, OU TPUT ,T APEi , TAPE2 ,TAPE3 ,T APE4 ,TA PE5 , 

1TAPE6= INPUT! 

DIMENSION  ARRAYC 4,1000 >,HV<3) 

PRINT  1 

1  FORMAT  <lHl,10X,MINE*,16X,*X*,24X,*Y*,24X,*Z*,////> 

DO  5  1=1,4 
00  5  J=i, 1030 
E  ARRAY ( I, J )=) • 

*  « 

*  WRITE  THE  HEADER  CN  TAPE5.  * 

*  ♦ 

READ  (fc,iO)  VR, DAY, SEC 
1C  FORMAT  < 3F12.0 ) 

WRITE  <5 )  YR , DAY , SEC 
20  N=u 

30  N=N*1  I  IFI LE=4 

*  * 

*  READ  THE  TRUE  DRAG  ACCELERATIONS  OFF  CF  TAPE4.  * 

*  * 

REAO  (O  <ARRAY(I,N>, 1=1,4) 

IF  (E0F(4)J  70,40 

•*  * 

^  REAO  THE  CORRESPONDING  THRUSTFR-I NDUCE D  ACCLERAT IONS  * 

*  ON  THE  X,  Y,  4N0  Z  AXES  OFF  OF  TAPES  1,  2,  AND  3,  RESPECTIVELY.  * 

*  * 

4o  DO  60  1=1,3 

IFI  LE=I  S  READ  (I)  T<II,HV(I) 

IF  (EOFmi  70,50 

so  j=m 


*  * 

*  TEST  TO  INSURE  THAT  THE  TIME  TAGS  MATCH.  * 

*  * 

TEST= A  RRA Y (l,N)-T (I) 

IF  (TEST. GE.. 06)  GO  TO  110 

*  » 

*  AOO  THE  THRUST-INOUCEO  VIBRATIONS  TO  THE  TRUE  ORAG  ACCELERATIONS.  * 

*  * 

6 j  ARRAY <J,N)=ARR AY <J,N)+HV(I) 

*  * 

4  ONCE  A  FULL  RECORD  OF  SIMULATED  DATA  HAS  BEEN  COMPUTED,  * 

*  WRITE  IT  ON  TAPE  5.  * 


IF  (N.NE.1JQ0)  GO  TO  30 

WRITE  (5)  N.ARRAV 

PRINT  eo,  ARRAY  <  GO  TO  20 

N=N-l 

IF  (N.EQ.O)  GO  TO  90 
WRITE  (5)  N, ARRAY 

PRINT  tO,  (CARRAYtl,  J)  ,1=1,4)  ,J*1,N, 

80  FORMAT  < l QX, F6. 0  ,5X, E2 0. 14, 5X, E20 . 14, 5X,E20 . 14) 
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GAUSS-MARKOV  PROGRAM 


PROGRAM  GAUSMK  C  IN  PUT, OUTPUT, TAPE4=INPUT ,TAPE1  ,T APE2 , T APE  3 > 

01  PE  NS  ION  A(17,4  > ,PM(3,3) , SM (3,3) ,PXI (3, 3) , RET At  3,3) , SETA  t  3,3) 

01  MENS  ION  WM(  3,1),WETA(3,1),WXI{3,1),W)3,1) ,PM1(3,1) ,PH2(3,1) 
DIMENSION  ETAOLO ( 3 , 1 > , PEI C 3, 1 ) , ET A  (3, 1 > , PE 2 T 3, 1) 

REAL  MM(3,i),  M(3,3) 

01  PE  NS  I  ON  PP1  (3,  1>  ,PP2  (3,1 ),  XI  13,  1»,XI0LDI3,  1)  ,BXI(3,  3) 

REAL  MP0L0(3,1) , MUB X , MUBY , MUBZ 

DIMENSION  SQ( 3,3>,PR (3,3) , RPR ( 3 , 3 ) , RNONOS ( 1 2) , ACC  1 100 C«4) 

REAL  I  SO (3, 3) , IRISH (3, 3) , IRESE(3,3> 

DIMENSION  RXI  (3,3) ,SW( 3,3) ,$(3,31,8(3,1) 

REAL  11(3,3), IRM(3,3),IRSE(3,3) 

DIMENSION  PRPRI3.3) ,CXII3,3) , ATRU (3 , 1 ) , ASEN ( 3, 1) , El ( 3, 1) , E 2 (3 , 1) 
DIMENSION  E3( 3,11 ,E4(3,1) 

PI=AC0S(-1.) 

7ERQ=0  . 

*  * 

*  WRITE  THE  HEADER  ON  TAPE  2.  * 

*  • 
REAO  (1)  YR,0 AT, SEC 

WRITE  (2)  YR, DAY, SEC 
PRINT  443 
443  FORMAT  (1H1) 

•  * 


♦  READ  THE  INPUT  CARDS. 

» 


* 


* 


DO  2C  1=1,17 

REA0(4,1C)  (A  (I,  .)>  ,J=1,4> 

10  FORMAT  (A10,E20.5,E2Q.5,E2C.5> 

PRINT  44b,  (A(I,  1)  ,  J=l,4) 

446  FORMAT  ( // ,5X , A1 9, E20. 5,E20.5 ,E20. 5) 

20  CONTINUE 
DO  33  1=1,3 
00  40  J= 1 , 3 

M(I,J)  =  RMI,J)=SM(I,J)  =  PXI  (I,J)=RETA(I,J)=SETA(I,J)  =  RXl(I,J)  =  0. 
SW(I, J)=0. 

40  CONTINUE 
30  CONTINUE 

II  (1, 1)  =  II(2,2)  =  II ( 3, 3)  =  1. 

II (1,2 1=11(1,3) =11(2, 1)=II( 2, 3) =11(3,1 )=IIC 3,2 )=0. 


*  ASSIGN  THE  INPUT  VALUES  TO  THE  APPROPRIATE  VARIABLES.  * 

•  * 

DO  50  1=1,17 

IF  (A(I,1).NE.10HBTASX  )  GO  TO  51 
HUBX=A(I,2)  8  VA RBX  =  A ( I , 3)  $  GO  TO  50 

51  IF  ( A ( 1 ,  1) .NE. 10HBIASY  )  GO  TO  52 

Muev=A(I,2>  3  VA  R8  Y= A ( I, 3)  $  GO  TO  50 

52  IF  (AU.n.NE.lflHBIASZ  )  GO  TO  53 
MU EZ=  A ( I , 2)  8  VARBZ  =  A(I,3)  $  GO  TO  50 

53  IF  ( A ( 1 , 1) ,NE • 10HT AUM  )  GO  TO  54 

TAUHX=  A ( I ,2)  $  T  AUMY= A  (1,3)  S  TAUMZ=A(I,4I  t  GO  TO  50 

54  IF  (A(I,1).NE. 10 HSIGMAM  )  GO  TO  55 

SM(1,1)=A(I,2)*PI/1A0.  t  SM(2,2)=A(I,3)*PI/180. 
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SH(3,3)=A(I,4)*PI/1«0.  S  GO  TO  SO 

55  IF  <A(I,  D.NE.IOHINTINT  1  GO  TO  5$ 

DELTAT=A(I,2>  S  GO  TO  50 

56  IF  (A<Itll.N£.lGHT AUET A  I  GO  TO  57 

TAUE  T  A  X=  A  { I,  ?  )  X  TAUFTAY=A  (I,  3)  l  T  AUE  T  A  7=A  (I  ,  4)  S  GO  TO  50 

57  IF  (AU.n.ME.inHSTGETA  )  GO  TO  58 

SETA  II  ,1>=A(T  ,2)  «  SET  A  (2,  2i  =  A  d,3>  i  SCT  A(3,3)=A(I,4)  5  GO  TO  50 

5fi  IF  (AT  I,  1>.NE.  IjMTAIJXI  >  GO  TO  59 

TAUXIX=A  (1,21  *  rAUXIY=A<I ,3)  5  TAUXI Z=A ( 1 , 4 )  t  GO  tq  50 

59  IF  <A<r,  11  .ME.  I0HPXI1  )  GO  TO  60 

PXIU,1)=ACI,2)  S  GO  TO  50 

60  IF  (A (  I, H .ME. 1GHPXI2  >  GO  TO  61 

PXl(i,2)=°XI(2,l)=Atl  ,2)  3  PXT ( 2 , 21=A ( I, 3 )  t  GO  TO  5D 

61  IF  <MI,U  »N".10HPXI3  )  GO  TO  62 

PXI(1,3)=PXI(3, 11=8(1, 2)  %  PX-I(3,2)  =  PXI  (?,3)=Att,3) 

PXK3,  3)  =8  (I,  4  )  5  GO  TO  50 

62  IF  <A(?,i>.NE.tOMSIGWAW  )  GO  TO  63 

SWli,l>=4Ct,?)  *  EW(2,2)=A(I«3)  S  S  W  <  3*  3  )  =  A  ( 1 , 4)  *  GO  TO  50 

63  IF  (A ( I, 1) .NE.10HSCAFAC1  J  GO  TO  64 

S(l*t)=A(I,2)  B  SI 1,2>=A(I ,3)  3  S 1 1, 3)  =  A ( 1 , 4  I  f  GO  TO  50 
6**  IF  (AIMI.NE.  10HSCAFAC2  »  GO  TO  65 

S  (2, 11  =  A  ( t,2)  «  S( 2*2)- All  ,3)  $  S(2,3)=A (1,4)  S  GO  TO  50 

65  IF  (A (I, 1).NE. 13HSCAFAC3  )  GO  TO  66 

S(3,1)=A(I,2>  3  S(3,?)=A<I,3>  S  S  T  3,  3>*  A(  I,  fc>  %  GO  TO  53 

66  IF  (A<r*l)*en, 13HTIMFLIM  )  GO  TO  67 
PRINT  <444,  A(I,1) 

444  FOfiiAT  ( 1H1,5X**FPR0P  ON  OATA  CAROS - IO  =  *,A1C) 

GO  TO  666 

67  TI"ELIM-A(I,2 ) 

53  coniNue 

*  * 

*  COMPUTE  THE  9 1 A?  FCP  El.  » 


CALL  T3M3FT(5fl,> 

CALL  TaM4MG (9  NON  OS , 3) 

911*1)  =  S'.,CT  (V  A*3  <)  *R  NO  NOS (  1)  +aU3X 
9(2,  1)  -  SORT (V  AR9  »)  *RNCNOS (?)  ♦'*U9T 
9(3,1) =  SCRT IV AR9Z) *RNONOS( 3) ♦  MU9Z 


« 

* 


OERIVE  SORT ( I -RH2) SH  F CR  E2.  * 

« 


RH  (l,11  =  FXP(-l.*OEL  TAT/TAUHX) 
PM(2,2)=CXP<-l.*0f LTAT/TAUMri 
RM (3, 3 )  =  FXP (~ 1 •*  OFLTAT /TAUM2) 

CALL  H PRC 0(3,  3«3,R**,3,RM,3,SQ,3) 

CALL  1SUPT(3, 3,II,3,SC,3,ISQ,3) 

OO  3J  1=1,3 
00  93  J= 1 » 3 

IRECI, J)=SQRT (ISOII, J)> 

93  COKINUE 
8o  COTTINUE 

CALL  HPRCO(3,3,3 ,IPM,3,SM, 3, IRMSM.3I 


*  DERIVE  RET A  FOR  ETA  OF  EX.  • 

*  * 

RETA(l,l)=EX<M-l.MOELTAT/T»UET»X>> 

RE  TA  ( 2 , 2  l=EXP  (-!,*( DEL  TAT/TAiJETAY)  > 
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RET  A (3, 3 )  =  EXP (  -1.*  (DELTAT/TAUETAZ) > 

*  DERIVE  SORT ( I-RETA2) SETA  FOR  ETA  OF  E3. 

* 

CALL  MPROO(3,3,T,RETA,3,RETA,3,SQ,  31 
CALL  MSU0T  (5,3, II , 3 , SQ , 3, ISO, 3) 

00  10  0  1=1,3 
00  110  J= 1 , 3 
IRSE (I , J)  =  SQRT (ISO  (I, J)) 

110  CONTINUE 
100  COtTINUE 

CALL  NCR 00 (3, 3, 3, IRSE, 3, SETA  ,3,IRESE,3> 

* 

*  DERIVE  RXI  FOR  XI  OF  E3. 

* 

RXI  ( 1,  1)=EXP(-1.  *0£LTAT/TAUXIX> 

RXI (2, 2)=FXP(-1. *DFLTAT/TAUXIY) 

RXI <3, 31 =EXP(-i. *0FLTAT/TAUXI7> 

CALL  MPR00(3,3,3 ,P  X 1 , 3  ,  RXI ,  3  ,  PP ,  3  > 

CALL  Mt5RC0(3,3,3  ,RX  I  ,  3 ,  PR,  3,  R  PR  ,  3) 

CALL  MSUBT (3,3,PXI , 3, RPR , 3 , PRPR, 3 > 

* 

*  FINO  3X1  ANO  CXI  USING  SQUARE  ROOT  ALGORITHM. 

* 

IF  (PKI(1,1I.NE.J.)  GC  TO  3H 
00  36  1=1,3 
DO  37  J=  1*  3 
CX  !(I, J>=BXI(I,JI=C, 

37  CONTINUE 
36  CONTINUE 
CO  TO  39 

33  CXI(1,2>=CXI( 1, 3>=CXI(2,3>=0. 

cxm,  u=sqrt  (pxi(iti)) 

CX 1(2,1 )=PXI( 1,2* /CXI (1,11 
CXI (3, 1)=PXI( 1,3) /CXI (1,11 
CX  1(2,21  =SORT  (PXH  2,2)  -CXI  ( 2  » l ) **2) 

CX  1(3,2)  =  (PXI (2, 3) -CXI <2,1)*CXI( 3,1) ) /CXI (2, 2) 

CXI (3, 3)=SQRT (PXI (3, 3) -CXI (3,1  >**2-CXl (3, 2) **2> 

9X1 (1,2) =8X1 (1,3 >=9X1(2, 31=0. 

9X111,1  )  =  SQRT  IPRORd,!)) 

9X1 (2, 1)=PRPR (1, 2) /BXI  (1, l) 

BXI (3, 1)=PRPR (1,3) /BXI (1, 1) 

3X1(2, 2)=SORT(PRPR(2,2)-9XI(2, 1)**2> 

3X I( 3, 2) =(PRPR (2  ,3) -9X1 (2,1) *3X1  (3,1)) /BXI (2, 2) 

9X1(3, 3 ) =SQRT (PR°R(3,3)-BXI(3,1)**2-BXI(3,2)**2) 

* 

*  PRINT  THE  INPUT  ANO  COMPUTFD  TIME-INVARIANT  MATRICES. 


* 

* 


* 

* 


* 


* 


* 


39  PRINT  901,  (3(1,1)  ,1=1,3) 

901  FOFMAT  ( 1H1 , 1  OX , *B I  AS  MATRIX* , /// , 10X, F2 0 . 5 , 1 0 X, E20. 5, 10 X, 
1E20.5) 

PRINT  902 

902  FORMAT  (//////////, 1C  X ,* S  MATRIX*,//) 

00  1001  1=1,3 

1001  PRINT  1000,  <S (I , J) , J=l,3) 

PRINT  903 

903  FORMAT  (//////////, 1 C X , *RM  MATRIX*,//) 


‘■'•W 


l 


t 


i  ' 
; 

& 

4) 


00  1002  1*1.3 

1002  FRINT  1000,  <RM< I , J) , J=l, 3) 

PRINT  904 

904  FORMAT  <//////////, 1QX ♦ *SM  MATRIX*.//) 

00  1303  1*1*3 

1C 03  FRINT  1C  CO.  <  SMt I ,  J) .3*1.3) 

PRINT  qC 5 

905  FQFNAT  (//////////, 1 GX ,* TR*SM  MATRIX*,//) 

DO  1004  1=1,3 

1C  G  4  PRINT  100C,  ( I°M0N(I, J), J=l, 31 
PRINT  9  05 

906  FORMAT  (//////////,  1GX,*RETA  MATRIX*,//) 

00  10 *5  1=1,3 

1005  PRINT  103C,  (RFTAd.J).  J*l,3) 

PRINT  907 

937  FORMAT  <//////////,  10X  ,  *SET  A  MATRIX*,//) 

DO  H36  1  =  1,3 

1006  FRINT  JOOC,  <  SET  A  (  I,  J)  .  J=  1 »  3) 

PRINT  PQ* 

«08  FORMAT  <//////////, 1GX,*IRFSE  MATRIX*,//) 

90  1007  1=1, T 

1007  PRINT  1030,  CIRESE(I.J».J=1,3> 

FRINT  703 

7C3  FORMAT  (/✓////////. 10X,*PXI  MATRIX*,//) 

DO  704  1=1,3 

7G4  fs>IMT  10CC,  JPXltl,  J>  ,  J,l,3) 

IF  (PXTtl.n.NE.  0.  )  GO  TO  450 
FRINT  451 

451  FOFMAT  <//////////> 

45E  fqImAT^ tlXt5H*****,5X,*STNCE  PXI(l.l)  EO'IALS  7ERO,  ALL  ELEMENTS 
1F*,5X,5H*****,/. 1X,5H*****,6X,*TME  BXI  ANO  CXI  MATRICES  MERE  SET 
UX,*TO  2EFO*,5X,5H*****) 

45 C  PRINT  P09 

909  format  (//////////, 10X, *9X1  MATRIX*,//) 

DO  100»  1=1,3 

UQ«  PRINT  lGwC.  <  BXI  (I ,  J)  ,  J=1 , 3) 

PRINT  701 

701  FOFMAT  t //////////, 1  OX , *CXI  MATRIX*,//) 

DO  130a  1=1,3 

lvC9  PRINT  103:,  tCXI <1 , J) ♦ J=l, 3) 

PRINT  732 

702  FO  IMAT  (//////////,  ItJX  ,  *SH  MATRIX*,//) 

00  1G1C  1  =  1,3 

1CU  PKIMT  lOOC,  <SH(I,  J),J=1,  3) 

100C  FORMAT  (/, 10X ,E2 3.5,1CX,E20.5, 1CX,F20.5) 

► 

»  READ  THE  INPUT  DRAG  ACCELERATION. 


* 

* 

* 


705  FORMAT  <1H1,1GX,49N*****  ALL  ACCELERATIONS  ARE  IN  IG)S 

1**) 

NOACEL  =  0 
PNCSET  =3 . 

llll-  READ  <11  N,  <<  ACC  <  J,  L>  ,L  =  1.4>  ,  J=1,N) 

IF  <EOF(D)  66b,  25 C 


### 
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*  j 


2 5„  00  200  ICNT=i,N 
TI=4CCtICNT,ll 
A T RU ( 1 « 1 )  =  ACC ( ICNT , 2) 
AT  F'J  (  2 , 1 1  =  ACC  ( IC  NT  ,31 
ATR'J(3,1I  =  ACC  (ICMT  ,4) 


*  COUNT  THE  NUMBER  OF  TIME  LIMES. 

NOACEL  =N0ACEL ♦ 1 

* 

*  CHECK  FOR  THE  TIME  LIMIT. 

TF  (TI.GT.TI-ElIMI  GO  TO  667 

* 

*  CONVERT  INPUT  ACCELERATIONS  TC  UNITS  OF  6S. 

* 

DO  13  1-1*3 

MRJ  ( I , 1  J  =  ATRU  (2  »  1 1  /  (  .0  098  0621) 

33  CONTINUE 

IF  C T I .GT.20. 1  GO  TO  5  0D 
DO  521  1=1, J 
4SENII,1>=ATRU<I,1I 
661  COHIN'JE 
GO  TO  6  C  2 

* 

*  OBTAIN  THE  ELEMENTS  OF  THE  RANDOM  VECTORS 

*  FROM  A  GAUSSIAN  SEQUENCE  OF  NORMAL  VARIATES. 

* 

EC  J  RN0  3ET  =P NOSET  *43  * 

CALL  TBMSEMRNOSET) 

CALL  T9MRNGUN0N0S,12I 
DO  25  1=1,3 
NM(I,ll=RN&NOS(I» 

WFT4<I,ll=RNONOS (1*31 
MX  I ( I  *  1 ) iRNON  0  3(1*6) 

M( I, II rRNDNOS (1*3) 

25  CONTINUE 

* 

*  COMPUTE  El 

*  THE  R A NO CM  91  AS  ANC  SCALE  F AC TOR-INOUCEQ  ACCELERATIONS 

* 

CALL  McRCO(3*  3*1 ,S , 3 , A TRU , 3, El , 3 > 

CALL  MA00(3,1,U, 3, El, 3, El, 31 

*  COMPUTE  MM  FOR  E2. 

* 

IF  (NOACEl.NF. 1)  GO  TO  730 
CALL  Moooon,  3,1  ,SN,3,WN,3tMM,31 

GO  TO  731 

730  CALL  MFRCD(3. 3 , 1 ,RM, 3 , MMCL D,  3,  PM1 , 3) 

CALL  MPRCD (3, 3,1 , IRMS" , 3 , K F,3 , P*2 , 3 1 
CALL  MftOD(3,l,PMt, 3,PM2,5, MM, 31 

731  DO  MM  1=1,3 
MMCLOCI,  II =MM ( I ,  II 

69  CONTI NUE 


COMPUTE  M  FOR  £2 


M(1,2)=MM(3,1» 
N(1,3»=-1.*MM(2,1I 
MI2, 1) =-t.*NM  13, 1J 
M  12, 3I  =  KH(1,1I 
Ml  3*  H =MM(2,1 ) 
M(3,2l=-l.*NMtl, 1) 


COMPUTE  E  2 

THE  TIMF-OEPENOENT  ORIENTATION  ERROR 

CALL  MPROOI3,  3,1  ,M,3,  A TRU, 3 , E2 , 3 1 

COMPUTE  ETA  a  NO  XI  FOR  E3. 

IF  (NOACEL.NE. 1)  GO  TO  T32 

CALL  MPRCDI3,  3,1, SETA, 3,  WET A , 3 , ET A , 3 > 

CALL  MPR CO  1 3, 3,1, CXI, 3, MX  I, 3, VI, 3) 

GO  TO  733 

732  CALL  MPRCO(3,3,l,RETA,3,ETAOLP  ,3, PEI, 3) 

CALL  MPRCO(3,3,l,lRESE,3,WFTA,3,P£2,3> 

CALL  MAOOI3,1,PE 1, 3 ,PE 2 , 3, ET A , 3> 

CALL  MPRC0I3, 3,1 ,RX I, 3 , XI 0 LO , ?, PP1 , 3 1 
CALL  MPRCCI3, 3, t, 8X1,3, MXI, 3, OF 2, 3 1 
CALL  MA00I3,1,PP1,3,PP2,3,XI,3> 

733  00  68  1*1,3 

ETAOLO  <I,1»=FTA< I,l> 

XICLO (1,11 =XI II, 1) 

68  COMINUE 


COMPUTE  E3 

THE  AOOITIVE  CORRELATED  NOISE 

CALL  MA00<3,1 ,ET  A, 3, XI, 3, £3, 3) 

COMPUTE  E4 
THE  WHITE  NOISE 

CALL  M°RCO (3,3,1 ,SW«3,w,3*E4,3) 

C0MoUTE  THE  RESULTANT  CORRUPTED  DRAG  ACCELERATION* 

CALL  MA00(3,1  ,E1 , 3, E2, 3 , ASFN , 3> 

CALL  N AOC (3,1 , E3 , 3 , ASEN, 3, ASEN, 3) 

CALL  MAOOI 3,1 , £4 ,3 , ASEN, 3, ASFN, 3 » 

CALL  MA00(3,1,ATRU,3, ASEN, 3, ASEN, 31 
502  PRINT  467,  NO  ACE L , I ATRU 1 1 , II , 1= 1 , 3 > 

467  FORMAT  I 1H0 , IX ,1 H* , 14 , 1H*,2CX,*TRUE  AX  =  *, E15.5 , 6X, *TRUE  AY  «  *, 
1E15.5.6X  ,»TRUE  AZ  =  *,E15.5I 

PRINT  46?,  TI,(ASEN(I,1I,I=1,3I 

468  FORMAT  (1X,*TIME  =  *,F15.5,4X,*SENSED  AX  *  *,E15.5,4X, 

1 ’SENSED  AY  *  *,E15, 5, 4X, ’SENSED  AZ  a  *,E15,5) 
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CONVERT  THE  ACCELERATIONS  TO  UNITS  OF  KM/SEC2  AND  STORE. 
00  54  1=1,3 

ASEN(I , 1>=ASEN(I,1)*(.CG980621> 

34  COFTINUE 

4CC(ICNT,2)=ASEN 11,11 
ACC(ICKT,3)=ASEN(2,H 
ACC<ICNT,4)=ASEN (3,1) 

20J  COKINUE 


WRITE  ONE  RFCORC  OF  SIMULATED  OATA  ON  TAPE  2. 

WRITE  (2)  N,« (ACC( J,L > ,L-1,4> , J=l, 1000) 

GO  TO  mi 

666  STOP  “EOF  ON  ACCELERATIONS  FILE" 

667  STOP  "TIME  LIMIT  REACHED" 

ENO 
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Defense  Mapping  Agency 
Aerospace  Center 
St.  Louis,  MO  63118 
ATTN:  Robert  Ballew 

Defense  Mapping  Agency 
U.S.  Naval  Observatory 
Bldg.  56 

Washington,  DC  20360 
ATTN:  Charles  Martin 

Defense  Mapping  Agency 
Topographic  Center 
Washington,  DC  20315 
ATTN:  R.  Smith 

Library  of  Congress 
Washington,  DC  20540 
ATTN:  Gift  and  Exchange  Divi¬ 
sion 
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Cameron  Station 
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